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Abstract — We propose a direct reconstruction algorithm for 
Computed Tomography, based on a local fusion of a few 
preliminary image estimates by means of a non-linear fusion rule. 
One such rule is based on a signal denoising technique which 
is spatially adaptive to the unknown local smoothness. Another, 
more powerful fusion rule, is based on a neural network trained 
off-line with a high-quality training set of images. Two types of 
linear reconstruction algorithms for the preliminary images are 
employed for two different reconstruction tasks. For an entire 
image reconstruction from full projection data, the proposed 
scheme uses a sequence of Filtered Back-Projection algorithms 
with a gradually growing cut-off frequency. To recover a Region 
Of Interest only from local projections, statistically-trained linear 
reconstruction algorithms are employed. Numerical experiments 
display the improvement in reconstruction quality when com- 
pared to linear reconstruction algorithms. 



I. Introduction 

THE FILTERED Back-Projection (EBP) algorithm is ex- 
tensively used for image reconstruction in the Computed 
Tomography (CT). This algorithm implements in the 2-D case 
a discretization of the inverse Radon transform. Despite the 
popularity of this method, its drawbacks are non-negligible: 
EBP fails to account for the numerous physical phenomena 
present in the data acquisition process (resulting in the well- 
known streak artifacts) and suffers from discretization errors. 
Moreover, it lacks the flexibility required to process partial 
input data, like in the case of projections truncated to a 
Region Of Interest (ROI), or in the case where projections 
are restricted to a limited angle range. 

Numerous techniques were developed to improve the per- 
formance of EBP. Many of them modify the filters, applied to 
the projection data. The basic problem with the standard Ram- 
Lak filter [1] is the high-frequency noise amplification. It is 
commonly treated by using an additional low-pass filter, which 
cut-off frequency is compatible with the expected bandwidth 
of the signal. Such an approach requires tuning the cut-off 
frequency and other parameters of the low-pass filter. The 
work reported in [2], for instance, is dedicated to tuning these 
parameters for lesion detectability. 

One important task in the clinical CT reconstruction is 
to recover the CT image in a Region Of Interest in the 
patient's body using low-exposure scan. To that end, there 
exist techniques which require projections data only in the 
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neighborhood of the ROI, in addition to a small number of 
full projections. These techniques employ an EBP modified to 
compute wavelet coefficients of the sinogram or the sought 
image by applying wavelet ramp filters [3]-[5]. This allows 
to recover the high spatial frequencies in the ROI from local 
data. The low frequencies require full projections; however, 
since they demand much smaller angular sampling rate, a local 
reconstruction can be done with much less X-ray exposure. 
There are, however, no proposed modifications of EBP that 
undertake the ROI reconstruction in the setup where only 
projections in the neighborhood of the ROI are measured. 
Such algorithm would allow further reduction in X-ray dosage, 
reduce scan duration and dismiss the necessity for the regis- 
tration of the projection data. 

To improve the reconstruction quality and overcome the 
above drawbacks, statistically-based iterative algorithms were 
developed. They use an elaborate model for the CT scan 
measurements, including sources of noise and partial projec- 
tions data. We refer to [6] for a detailed example of such an 
algorithm and to [7] for a broad overview of iterative methods. 
Unfortunately, high computational cost of these algorithms 
restricts their use in clinical CT scanners. 

Superiority of the statistically-based methods stems from 
the fact the image is reconstructed in a non-linear, locally- 
adaptive manner, relying only on the available projections 
data. Eor instance, with the method described in [6] the 
reconstruction is performed via optimization of a Penalized 
Weighted Least Squares (PWLS), with a penalty promoting 
local smoothness in the image domain in an edge-preserving 
manner. Such behavior can not be achieved with a linear, 
spatially-invariant algorithm like EBP. Therefore, an algorithm 
combining the advantages of both approaches is desired: a 
direct (and, therefore, fast) processing of the available data, 
on one hand, with a non-linear, locally-adaptive nature on the 
other hand. 

One such algorithm is developed in [8]. The proposed 
algorithm uses a powerful filtering technique, involving a 
training procedure. It employs an exemplar-based classification 
of the sinogram data patches, combined with training of local 
2-D projection filters, which results in a non-linear overall 
filtering procedure. 

In this work we propose a non-linear, locally-adaptive 
reconstruction scheme, based on example-based statistical 
learning of its components. The scheme consists of two stages: 
first, a sequence of linear EBP-like transforms are applied 
to the available projections data, resulting in a number of 
preliminary image estimates. Then, a local non-linear fusion 
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of these estimates is performed to produce the final image. 

In the setup of an entire image reconstruction from a full- 
scan, the linear estimates are obtained with the FBP algorithm 
with a varying cut-off frequency of the projections filter. A 
more important goal is to reconstruct an ROI from truncated 
projections. Since FBP does not perform well in the absence 
of global projections, we have developed a more flexible linear 
reconstruction scheme called AFBP [9] . It generalizes the FBP 
and employs more powerful filters in the sinogram and image 
domains. The convolution kernels for these filters are derived 
via statistical training which accounts for missing data and the 
desired reconstruction properties. In the proposed algorithm, 
a set of linear estimates of the ROI image are computed by a 
number of AFBP versions only from projections through a disk 
containing the ROI, which radius is 110% of the ROI radius. 
Then the ROI is recovered by a neural network from these 
preliminary reconstructions by a non-linear, learned fusion 
rule. 

While the computational cost of this procedure is only 
about ten times the cost of the FBP algorithrrQ, two main 
features distinguish the proposed method from any spatially- 
invariant reconstruction transform. First, the components of 
our scheme are designed to work with partial data (truncated 
projections). Second, the reconstruction is locally adaptive in 
the image domain, which allows to reduce the noise present in 
linear estimates and to preserve edges and texture in a better 
way. The proposed method is labeled as SPADES (SPatially 
ADaptive Estimator). 

This paper extends our work presented earlier in conference 
publications [9], [10]. Most of the presented material is new, 
except for the description of the linear AFBP scheme. The 
SPADES algorithm, which is also briefly presented in [10], 
has been revised and improved. 

The paper begins with Section [III containing preliminaries 
and notation. For a theoretical motivation of SPADES we 
present the locally-adaptive denoising algorithm of Lepski, 
Goldenschluger and Nemirovsky [11] and extend it to CT 
reconstruction setup (Section [IIl|). Then SPADES is described 
and demonstrated in Section HVl The AFBP scheme is intro- 
duced in Section |Vl Then, the SPADES is extended to ROI 
setup in Section IVll Discussion follows in Section IVlIl 

II. A MODEL OF 2-D Transmission Tomography 

A 2-D slice of a physical object is represented by the 
attenuation map f{x), defined in the domain A C R'^ - 
this is the image recovered by the CT reconstruction. This 
map, assumed to have a support radius R, is projected along 
straight lines by means of the 2-D Radon transform: for 
s e [—R^R],0 G [0,7r], the transform g = R/ is defined 
by 

9e{s) = / f{s-cos{0) — t-sin{0),s-sin{0)-\-t-cos{0))dt. 

JteR 

(1) 

^In practice, the run time can be reduced to one FBP computation by 
parallel execution on small number of cores. 



We denote the range of R by 7^ C M^. The adjoint transform 
R*, also known as a Back-Projection, is defined by 

(R^g){x) = f g0{[cos{e),sin{O)]-x)dO. (2) 
Je 

In the discrete setting, the Radon transform of f{x) is sam- 
pled at a large number of fixed angles (evenly covering the 
range [0,7r]) and fixed signed distances s (bins). The matrix 
image f{x) is computed from these samples by a discrete 
reconstruction algorithm. 

Let i = £{0^ s) be a line which makes the angle 7t/2 — with 
the X axis and passes at a distance s from the origin. To each 
such line there corresponds a detector which counts the num- 
ber y£ of photons in a specified time interval during the scan. 
Due to the limited photon count, the values y£ are modeled as 
realizations of random variables: m ~ Poisson (/oe"'^^-^^^). 
The X-ray source intensity Iq determines the parameters 
of the Poisson distribution, thus controlling the noise level. 
The Maximum Likelihood (ML) estimate of R/ from the 
measurements m is ge{s) = —log{y£/Io). In the ideal case 
where yi = ^(1^) = /oe~*^^-^^^ (not attainable in reality 
since the expectation of Yi is not, in general, an integer), ML 
estimate is the true Radon transform of f{x). 

More realistic modeling of the CT scan takes into account 
additional disruptions, such as electronic noise (additive con- 
stant in the Poisson parameter A^), scatter of X-rays (Compton 
effect), crosstalk among the detectors, and more [12]. 

The basic FBP algorithm is defined by means of the linear 
reconstruction operator Tfbp = FRam-Lak- Here the 

filter FRam-Lak IS the 1-D convolution kernel, applied to each 
projection ( acolumn in the sinogram matrix). The Ram-Lak 
kernel ( is defined in the Fourier domain by ({uj) = \uj\. As 
mentioned earlier, the Ram-Lak kernel is often smoothed by 
some low-pass filter compatible to the reconstructed images 
and the noise level. 

III. From Adaptive Denoising Technique to a 
Fusion Rule for CT Reconstruction 

Our goal is to bridge the gap between the linear and the 
iterative algorithms by a direct reconstruction scheme, locally 
adaptive to the data. To that end, we employ the technique of 
filtering with adaptive kernels, developed originally for signal 
reconstruction from noisy measurements. In the classical setup 
of the problem, a signal / is measured through 

where the set {xi} is a sampling grid and are independent 
normal random variables. The task is to compute an estimate 
f{x) which minimizes the L2 norm of the error e = \ f — f\ 
in an interval of interest. 

One basic denoising technique is to apply a linear kernel 
estimator = y ^ where /t: is a constant convolution 
kernel representing a low-pass filter. Such an estimator will 
mistreat those regions where the bandwidth of the signal does 
not match the bandwidth of the filter (i.e., either the edges 
will be blurred or the smooth regions will remain noisy). 

A substantial improvement can be achieved by using, at 
each image location, a low-pass filter which cut-off frequency 
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matches the local spatial smoothness of the signal. Thus, in 
smooth regions a stronger blur will be used, averaging the 
noisy values, and near edges almost no blur will be applied 
in order not to smear them. Of course, such knowledge of 
the local smoothness of the underlying clean signal is not 
available, but it can be evaluated using local image statistics. 
Therefore, given a sequence of filter kernels with gradually 
growing measure of the blur and a good decision rule, choos- 
ing an appropriate filter for each spatial location, a spatially 
adaptive signal estimator can be implemented. 

An analogue of image filtering with kernel Hi in CT re- 
construction is a linear reconstruction transform T with the 
property that the Point Spread Function (PSF) of the (shift- 
invariant) operator TR is the kernel hz. It is easy to obtain 
a sequence T^,z = 1,...,/ of linear reconstructors with 
gradually growing spread of the PSF; therefore an analogue 
of spatially-adaptive denoising can be designed for CT re- 
construction, as we detail in Section IIII-BI We should stress 
that the algorithm, described and demonstrated therein, is 
not intended for practical use: it is a mid- step between the 
theoretically-based denoising technique for Gaussian noise and 
the non-linear CT reconstruction algorithm based on a learning 
machine. 

A. Lepski-Goldenshluger-Nemirovsky Estimator 

A locally adaptive estimator for signal denoising was de- 
vised in the work [11] by A. Goldenshluger and A. Ne- 
mirovsky, based on a general scheme by Lepski. Consider 
again the task of recovering a 1-D signal f{x) from noisy 
observations y{x) = f{x)-\- ^(x) on a discrete grid points 
{x G r}, where {£,{x)}xer is a sequence of independent 
normal random variables. To estimate f{x) at some point 
xo, a linear combination of the neighboring samples y{x) is 
taken (as is done by a convolution kernel with corresponding 
coefficients). The basic idea of this estimator is demonstrated 
for simple rectangular windows centered in xq'. f{xo) = 
mean{y{x) \x G A}, where A = [xo — S^xo-\-S]. In practice, 
order-m Least Squares (LS) polynomial is fitted locally into 
the data. 

The error |/a(^o) — /(^o)| consists of the deterministic 
dynamic error ujf{xo, A) = meanxeA\f{x) — f{xo)\ and the 
stochastic error C(^) = ^XIxga^^(^)' introduced by the 
noise. As the window A grows, the dynamic error increases 
(with the rate related to the signal smoothness) and the 
stochastic error decreases, since the noise is averaged over a 
larger interval. The goal is to find an optimal window width, 
which balances the error components. 

In order to approximate the optimal width, the authors 
of [11] employ confidence intervals. Consider a sequence 
Ai, Aat of windows, centered in xq, with a growing width. 
For each A^ there is an estimated upper bound pi of the 
stochastic error in this window. For i = 1 it is assumed the the 
window is small enough so that the dynamic error is dominated 
by the stochastic one, hence the overall error is bounded from 
above by 2pi. The confidence interval related to A^ is then 
defined by 

A = [/Ai(^o) - 2pi,/Ai(xo) + 2pi]. (3) 



Notice that for every index i satisfying cj/(xo, A^) < Ci^i)^ 
the total error |/Ai(^o) — /(^o)| is smaller than 2pi (with the 
aforementioned high probability). Therefore, the interval Di 
contains the true value f{xo). 

In order to estimate the maximal index z* with this prop- 
erty, the intersection of intervals Di is considered. A simple 
argument, presented in [11], shows that if i+ is the maximal 
index for which 

n Ay^0, (4) 

then the estimate error is bounded by 

\fi+{xo)- fixo)\<Gpi' (5) 

Since the stochastic error pi* is minimal among possible errors 
corresponding to different segments A^, the choice of fi+ is 
nearly optimal. 

As the simple averaging is replaced with order-m LS 
approximation, the efficiency becomes comparable to best 
denoising techniques. It is proven in [11] that such algorithm 
is near-optimal. 

Before passing on to the CT reconstruction, we mention 
that a similar algorithm was devised by Lepski, Mammen, 
and Spokoiny in [13]. Both works implement a general 
scheme of Lepski [14] and seem to share the same approach. 



B. A Switch Rule for CT Reconstruction 

The denoising technique, described above, requires a pre- 
liminary sequence of linear signal estimates, obtained by 
filtering the noisy signal with a corresponding sequence of 
convolution kernels (rectangular windows of growing radius, 
in this case). In the setup of CT, we use a sequence of linear 
transforms {T^j^^^, each is the FBP algorithm involving a 
low-pass filter with a different cut-off frequency. Let g denote 
the noisy sinogram and fi = Ti{g) is the sequence of image 
estimates. The the noise, present in the reconstructed images, 
is of a complicated, data-dependent nature. The denoising 
algorithm requires to know the bound Pi{p) on the statistical 
error in each image location p of each these bounds can 
be computed using noise variance in the image domain. 

Estimating the noise variance is difficult, and we do not 
pursue this task here. Our goal is to show that given the 
necessary information for the switch rule, locally adaptive CT 
reconstruction is possible. Once the evidence for such success 
is obtained, we propose a different fusion rule, based on a 
learning machine. Therefore, we simulate a large number on 
noise instances to compute the variance \i{p) of the noise 
at the location p in the image estimate fi. Then, instead of 
estimating the bound pi on the stochastic error in p as it is 
done in [11], we compute 

Pi{v) = <VMp)Y (6) 

where the parameters i<i.,q are tuned using a grid search on 
the training set. These two parameters are required to calibrate 
the numerical values of the data obtained in CT reconstruction 
process and are set exactly once for the given setup. 



The confidence intervals are computed by the formula O 
Then a switch rule, expressed by the condition is applied 
to compute the index The output value in location p is 
fi+{p). We label the described algorithm as LeGoNe after 
the authors of the prototypical denoising technique (Lepskii, 
Goldenshluger, Nemirovski). 

C. Numerical Experiment - CT reconstruction with LeGoNe 

The quality measure we use is the Signal-to-Noise Ratio 
(SNR), defined for the true signal /o and its estimate / by 

SNR{foJ~) = -20log J^'~/^^\ (7) 

II/0II2 

In this experiment we use a set of randomly generated 
256 X 256 geometric phantoms. Each phantom constitute of 
a large ellipse with boundary and a constant background, 
filled with many smaller ellipses with randomly chosen centers 
and radii. There are four intensity levels in each image 
(also randomly chosen for each phantom). Examples of the 
phantoms are presented in Figure [T] 




Fig. 1: Examples of geometrical phantoms used in the exper- 
iments. 

In our simulations of CT projection and reconstruction, 
the projection sets (sinograms) of the reference images are 
computed using a pixel-driven implementation of a discrete 
Radon transform R. Its adjoint transform is employed in the 
reconstruction process. The sinogram noise is generated in 
accordance to the statistical model, described in Section III 
(see also the Appendix section for the details). 

The sequence of the EBP algorithms {Sijf^i is generated 
by applying a Butterworth window with growing radius to the 
Ram-Lak convolution kernel k,o in the Eourier domain: 

f^, = f^o^F-\H,), H = H{p,q,) = -^-^ (8) 

1 + q/ 

Here F is the 1-D Discrete Eourier Transform. The parameter 
p controls the steepness of the window roll-off, and q (cut-off 
frequency) determines its width. is a monotonously 

decreasing sequence. Figure [2] displays (part of) the corre- 
sponding sequence of fi = S^^ obtained from a sinogram g 
of a geometric phantom. 

The LeGoNe reconstruction algorithm, as described earlier, 
is applied using the images Comparison of various 

reconstruction results can be observed in Figure [51 Notice 
that, despite the modest increase in the SNR value, much of 
the noise present in EBP reconstructions is removed in the 
LeGoNe estimate. 

In the switch map generated by the LeGoNe (Figure 131), an 
intensity value in each location p is the index z+ chosen for the 
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Fig. 2: Sequence of linear EBP estimates for the projected 
phantom, with a growing degree of blur. 




Fig. 3: Reconstruction results, left to right: EBP with ideal 
Ram-Lak filter, EBP with optimally apodized Ram-Lak filter 
(SNR = 15.72 dB), LeGoNe (SNR = 16.13 dB). 



outcome image in p. It can be observed that in smooth areas 
of the phantoms the algorithm prefers higher indices (since 
in blurred images the stochastic error is lower) and near the 
edges of the ellipses, lower indices are chosen. 




Fig. 4: The switch map of the LeGoNe algorithm. 

In practice, the LeGoNe algorithm is not our method of 
choice. First, the underlying denoising technique is optimal 
only in the mini-max sense and up to a constant. Second, 
the evaluation of noise statistics in the image domain is 
difficult. More importantly, we can not use outputs of other 
reconstruction algorithms as preliminary results, to be further 
improved by LeGoNe. 

IV. SPADES - Local Fusion Based on a Neural 
Network 

To build a more powerful local fusion rule we resort to a 
neural network. Indeed, it is difficult to devise an analytical 
rule to approximate the true image value from the set of pre- 
liminary reconstructions fi = T^g. The problem is severed by 
the non-homogeneous, data- and algorithm-dependent noise, 
present in the image. Yet, by the analogy with the situation 
discussed in the previous section, we can hope that the set 
of values {vi = fiix)}^^^, for carefully chosen estimates 
fi, contains information on f{x) which is more accurate 
than any of the individual estimates; in other words, there 
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exists a function ^ = ^(^i, v/) that will produce a more 
accurate image reconstruction when applied point- wise to the 
set {h}U- 

A neural network can learn this function via a training set 
of reference images {f^}J^i and their simulated estimates 
{fi}i<i<i, i<t<T- The sinogram data is obtained by either 
applying a Radon transform to discrete reference images, or 
(in a more realistic setup) by scanning a set of geometric 
phantoms in a clinical CT scanner. Then the linear reconstruc- 
tion algorithms are applied in order to compute the training 
data {//}. The learning procedure occurs off-line, and when 
an unknown object is scanned, the images {fp}i=i are 
computed from its projections data and fed to the trained 
neural network point-wise, to produce the final outcome. 

Notice that in this setup we are not confined to use a 
sequence of linear reconstruction algorithms with a gradually 
increasing measure of blur, like with the LeGoNe technique. 
Any available preliminary reconstructions (not necessarily 
linear) can be fused with the neural network in order to 
produce a result, expectedly superior to all the participating 
versions. The scheme of the resulting SPADES algorithm is 
given in the Figure [5]). 

Linear Reconstruction 
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Neural 
Networl< 

/Feature \ ^^^1 f 

\Extraction/ ^^^H 




Fig. 5: SPADES reconstruction scheme 



A. Algorithm Description 

Definition of the neural network: We use a single-layer 
feedforward type neural network. Its output function is defined 
by 



N 



K 



(9) 



k=l 



where N is the number of neurons, K is the size of the input 
vector of features x, Wkj is the weight on edge connecting k- 
th input to j-th neuron and a{z) = z/{l-\-\z\) is the sigmoid 
function. 

The learning procedure consists in solving the following 
optimization problem: given a set of vectors of features = 
[x\^...^x\^] and the corresponding true values v{t) for t = 
1, ...,T, minimize the objective function 

T 

(k;*, 'U*) = arg mm{S^{v{t) — y{x^ ^w^ ^))'^} (10) 



This function is non-convex. In our experiments it is min- 
imized using the Matlab routine fminunc.m. It performs an 
unconstrained optimization using the BEGS Quasi-Newton 
method with a cubic line search procedure. 



Training set design: Notice that intensity values of the 
reference image f{p) can vary significantly over different 
regions. We reduce the variability of the data fed to the neural 
network by using relative values of the images: the vector of 
features Xp corresponding to a spatial location p is build as 
follows. The first / values are set to 

xA^=Up)-f{p). (11) 

where / is the besH available estimate of f{x). Other entries 
of Xp consist of image values in a small neighborhood J\fp of p 
taken from /. This provides the neural network with additional 
local information about the image at the point x, allowing for 
more accurate restoration of the value f{x). 

The corresponding true output value y^, provided along with 
the vector Xp in the training stage, is set to yp = f{p) — f{p)- 
In the reconstruction stage, the final outcome image is obtained 

fnn{p) = y{Xp, V, W) + f{p). (12) 

The dynamic range of the input values for the neural 
network is normalized to [0, 1]. 

B. Numerical Experiment - SPADES on a full-scan data 

We repeat the previous experiment in CT reconstruction, 
when instead of the LeGoNe switch rule the fusion is per- 
formed by a neural network. The training data is extracted 
in the way described earlier in this section. A network of 24 
neurons was trained with a set of 15800 vectors of features 
(sampled from ten images) and then applied to a test phantom. 
Ten preliminary versions of each image were built, using an 
EBP algorithm with varying degree of blur. 

The outputs of the EBP and the SPADES algorithms are 
displayed in Figure [6l It can be observed that the noise streaks, 
characteristic for EBP reconstruction, do not appear in the 
SPADES output. It is much closer to the piecewise constant 
phantom, and the SNR value reflects the quality improvement: 
it is 2.3 dB higher than the linear reconstruction and 1.8 dB 
higher than the LeGoNe reconstruction. 




Fig. 6: Left to right: phantom, best linear reconstruction (SNR 
= 15.72 dB), neural network output (SNR = 18.00 dB) 



V. AFBP Reconstruction Algorithm 

We describe a design of a linear reconstruction transform, 
introduced earlier in the conference publication [10]. 

^In the SNR sense. We use the training set to determine which linear 
estimator is the best. 
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A. Definition of the AFBP Transform 

The extended linear reconstruction operator, labeled as 
Adaptive Filtered Back-Projection (AFBP), is defined in the 
block diagram given in Figure [71 with parameter set n = 



Sinogram 
filter 


Back- 
Projection 


Image 
filter 






R* 











Fig. 7: AFBP reconstruction scheme. 



The linear filter Y^v acts on a sinogram by a distanceH- 
dependent 2-D convolution kernel, which definition is detailed 
below. The filter F^^ in image domain is implemented as a 
single 2-D convolution kernel applied to the output of the 
Back-projection transform.The involved kernels are generated 
in the training process, described in the Sectior TV-Bl 

Filter - using 2-D convolution kernels: The ideal 
inverse Radon operator only requires a 1-D filter applied 
to each projection. In the practical setup, we use a two- 
dimensional kernel applied to the sinogram (each projection 
is also affected by the few neighboring ones). Such a filter 
exploits the correlation between the neighboring projections 
and can improve the reconstruction quality. 

Filter Y^v - using a distance-dependent kernel: When 
the projections are truncated to the ROI, the central part of 
each remaining projection should be filtered using a symmetric 
kernel with a small spatial support in order to reduce the 
truncation error. Near the edges the information should be 
gathered in a non- symmetric way, only from the non-truncated 
part of the projection. This can be done by assigning a separate 
kernel to segments of radial distance from the center of the 
projection. Explicitly, assume that only projections available 
are over lines passing through a central disk of radius D in the 
image domain. We partition the range of distance s G [O^D] 
into d disjoint sub- segments [O^D] = UiLi ^se a bank 

of K,^ = {K,f}f^^ of corresponding convolution kernels. The 
filter is applied by convolving the two projection segments 
corresponding to ±Di with the kernel hzf . See Figure [8] for 
an illustration. 

Here we use the AFBP algorithm in the setup of partial 
data - truncated projections. This technique also improves FBP 
performance in the full- scan setup, but in this case the gain is 
not substantial. 

B. AFBP Parameters Training 

We state two goals pursued in the parameters training of 
the AFBP: 

Goal 1: Assume availability of a representative set from 
some family of images. X-ray intensity of the source, radius of 
a central disk where the projections are measured, and a radius 
of a central disk where the image should be reconstructed 
(ROI). The goal is to maximize the reconstruction quality of 
T, in the Mean Square Error (MSE) sense, for these images 

^we refer here to the distance s of a projection bin from the center of the 
projection 



Filtered by kappa-3 




Filtered by 
kappa-2 



Fig. 8: Illustration of the radially-variant sinogram filtering. In 
the left part of the diagram: concentric disks in image domain 
corresponding to horizontal bands in the sinogram. 

(the reconstruction is performed from noisy and truncated 
projections). 

Goal 2: Under the same conditions but regardless the X- 
ray intensity, design a reconstruction operator T with such that 
the action of the projection-reconstruction operator TR will 
be as close as possible to the action of a radially- symmetric 
Gaussian convolution kernel in image domain. 

The second goal follows from the needs of the non-linear 
fusion algorithm, which admits at its input a sequence of 
images reconstructed in different ways. A special case of such 
sequence, motivated by the LeGoNe denoising technique, is 
a sequence of images with a growing measure of blur. The 
AFBP training, corresponding to the stated second goal, will 
provide us with corresponding sequence of linear operators. 

Objective function 1: best image quality. The objective 
function addressing the first goal is designed as follows. 
Denote by r the ROI radius. Let Yrqi be the operator which 
nullifies all the bins in the sinogram matrix which correspond 
to line integrals outside the central disk (in image domain) 
of radius l.lr. Also, we denote by Froi an operator in the 
image domain that nullifies all the image values outside the 
ROI. 

Given a set Ttr of representative images, we build the 
training set Qtr consisting of noisy truncated sinograms 

Gtr = {g} = YROliKf + ef) I / e ^^r, j = 1... J}, (13) 

where {gj}j^i are obtained from / G Ttr by applying the 
Radon transform and generating J instances of Poisson noise. 
Using this training set, the parameter set k, of the AFBP 
transform is then computed as an optimizer of the objective 
function 

/.* = argmm ^ ||F^o/(T,^'/ - /) II2, (14) 
feTtrJ=i...J 

where is defined in previous Subsection. 

The objective ([T4l) is quadratic in parameter sets hz^ and in 
hz'^ separately. Thus, when fixing one, the other is updated by 
solving a corresponding linear problem (we use the Conjugate 
Gradients (CG) method). The training is then carried out in 
turns, fixing one set of variables and updating the other. We do 
not provide a convergence proof, but the property of the CG 
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algorithm guarantees a monotonous decrease of the objective 
function. In practice, we continue the update process until the 
first 5 significant digits of the objective function are stabilized. 
In our experience, no noticeable change in the algorithm 
behavior occurs in case of further optimization. 

The resulting linear reconstruction transform T,^ has the 
following advantages over the conventional FBP: 

(1) The objective function adapts the reconstruction operator 
to partial data conditions, and thanks to its distance-dependent 
filter such adaptation results in more adequate treatment of the 
truncated projections. 

(2) The adjustment of the low-pas filter parameters in FBP 
is replaced with automatic derivation of the filter, which cut- 
off frequency adjusts to the typical spectra of training images 
and the noise intensity. 

(3) The post-processing filter F^a is matched to the pro- 
jections filter. It helps reducing the value of objective function 
substantially (as compared to using only projections filter) and 
therefore to improve the output quality. 

(4) The specialization of the reconstruction operator to the 
given set of images makes it less universal, but improves the 
reconstruction of similar images. By using a set of images 
typical to the given task (for instance, sections of CT scans of 
specific body part), a dedicated reconstruction transforms can 
be trained for different clinical cases. Another application is 
to devote a personal reconstructor to an individual patient and 
use her previous scan data (say, from older similar CT scans) 
to allow lower X-ray dosage in future scans. 

Objective function 2: approximate convolution filter. The 

objective function built for the second goal involves noiseless 
projection sets. It has one additional parameter - the standard 
deviation a of the Gaussian kernel which action TR should 
mimic. The parameter set k, corresponding to the value of a 
is computed as the optimizer of the following equation: 

/.* = argmin V ||(T,R/ - F^o/(G, * /))||i (15) 
feJ^tr 

This objective function encourages the reconstruction trans- 
form T to produce a blurred output image, which resembles 
as much as possible the "true" scanned image, filtered with the 
corresponding Gaussian. Such reconstruction is an alternative 
to using an FBP algorithm with a low-pass filter applied on 
the projections. Advantage of the AFBP transform consists 
in the points mentioned earlier and allows us to use the non- 
linear fusion in the setup of ROI reconstruction from truncated 
projections. 

Trained kernels: We display visual examples of traned 
convolution kernels for the sinogram filter of AFBP. in Figure 
[9] displays two instances of such kernels, one trained for sharp 
PSF of the TR operator, and the other one for a wide- spread 
PSF. The central column of each kernel resemble the Ram-Lak 
filter, and the neighbor columns ( corresponding to neighbor 
projections) are similar except for the central main peak.The 
bin values of kernel producing sharp PSF decay rapidly away 
from the center of the kernel, and the bin values of the blurring 
kernel are non-decaying, providing input from far projection 
bins in the convolution process. 




(a) Kernel for sharp PSF (b) Kernel for wide-spread PSF 







i 


i 



(c) Ram-Lak kernel 

Fig. 9: Bin values of the trained kernels, compared to the 1-D 
Ram-Lak kernel. 

C. Comparison to Previous Work 

We return to the previously mentioned work reported in [8]. 
The Nonlinear Back-Projection (NBP) algorithm, proposed by 
its authors, employs a different kind of a non-linear filter. It 
is applied locally to the sinogram data, pre-filtered with the 
apodized Ram-Lak filter, in the FBP framework. Each small 2- 
D sinogram patch is filtered by a data-dependent combination 
of pre-learned filters. The parameters of these filters are sta- 
tistically trained to yield the lowest MSE values in the image 
domain for the training set of images. The proposed approach 
elegantly employs the Gaussian Mixture model in the sinogram 
domain and its patch-wise classification of the sinogram data 
has a very good potential to incorporate the sinogram features 
that would allow high-quality reconstruction. 

This technique shares it conceptual approach with the 
AFBP and the SPADES algorithms, proposed in this article, 
but can not be directly compared to either one. AFBP also 
employs a spatially-variant, statistically trained 2-D filter in 
the sinogram domain, but this filter is linear and not data 
dependent. SPADES is more closely related to NBP since 
both algorithms are locally adaptive: NBP employs a sort 
of data-dependent fusion in the sinogram space, by means 
of a combination of local filters, and SPADES acts in the 
image domain, fusing the pre-computed linearly reconstructed 
images with a neural network. It is difficult to argue in 
favor of one or the other approach, except maybe for the 
fact SPADES uses weaker assumptions on the reconstruction 
problem (we do not assume Gaussian mixture model on small 
sinogram windows). Unfortunately, we weren't able to conduct 
a numerical comparison of the two techniques: the results on 
geometrical phantoms, presented in [8], are inferior to FBP, 
and the successful experiment on the clinical SPECT data can 
not be repeated for technical reasons. 
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D. Numerical Experiments: AFBP on ROI Data 

We present numerical experiments of ROI reconstruction 
from truncated projections. The first one is performed with 
geometric phantoms, described in Section [Till The projections 
of 15 training images were truncated to a central disk of radius 
32 + 3 pixels. A Poisson noise, expressed in restriction of the 
source intensity to /q = 1200, was applied to the data. The 
ROI consists of a central disk of radius 32 pixels. The AFBP 
projections filter is a 2-D radially- variant convolution kernel, 
as described in Section IV-Al Specifically, the size of each 2-D 
kernel is 5 x 72 [angles x projection bins], and the ROI radius 
is divided into 5 sub-segments Di (overall, 5 2-D kernels). The 
parameters of the AFBP algorithm were trained using equation 
([T4l) and then applied to the truncated, noisy projections of the 
test image. 

For comparison we use the FBP algorithm which is imple- 
mented as follows: (1) The Ram-Lak kernel is smoothed by 
applying a Butterworth apodization window in the frequency 
domain. The parameters of the window function are tuned 
on the training set to provide the best SNR value. (2) The 
truncated projections are extrapolated by replicating the first 
and the last non-zero rays (see Figure [15]). This simple but 
effective linear sinogram completion technique was described 
in [5]. See Figures IIQIIII for a visual comparison of the 
algorithms. 




Fig. 10: Reconstruction in a ROI. Upper left: true image, 
upper right: AFBP reconstruction from truncated projections 
(SNR = 18.44 dB). Lower left: FBP reconstruction from 
truncated projections (SNR = 14.55 dB). Lower right: FBP 
reconstruction from full data (SNR = 18.06 dB). 

The experiment was repeated on a set of clinical CT images, 
which represent axial head sections (courtesy of the Visible 
Human project) The CT scan was performed on a GE scanner 
and consists originally of 512 x 512 grayscale images with 
pixel depth of 12 bits. The images were cropped and resized 
to 256 X 256 size for the purpose of our experiments. Some 
of them are displayed at Figure [T2|. Figure [13] displays the 
ROI and a typical result of reconstruction from truncated 

4 www. nlm . nih . gov/research/visible/visible_ human . html 




Fig. 11: Error images \forig — /| in the ROI. Left to right: 
AFBP on truncated projections, FBP on truncated projections, 
FBP on full data. Darker shades correspond to larger error. 



projections. In Figure [14] we compare the AFBP versus the 
FBP algorithm. 




Fig. 12: Clinical CT images from the AFBPt training set 



Fig. 13: Left: true image with a marked Region Of Interest. 
Right: full reconstruction from truncated projections by AFBP. 
The ROI is clearly distinguished. 

It can be seen from the images that AFBP improves sub- 
stantially the FBP output, producing an image comparable to 
a reconstruction from full scan data. The SNR values support 
this observation. 



E. Study of AFBP and FBP Algorithms on ROI 

We now turn to second kind of AFBP transforms T^, trained 
by optimizing the equation ([T5]) for different values of a^. 
These reconstruction transforms come to replace a sequence 
of FBP transforms {S^jf^^ in ROI reconstruction, due to 
the need to adjust the linear reconstructors to the missing 
data setup: FBP is theoretically derived for full-data scan and 
its performance of truncated projections is poor, even with 
optimized parameters of the low-pass windo\\|j. 

Now we define a measure on a reconstruction transform 
X that expresses the amount of blur X introduces into the 
output image. Let Go- stand for the rotationally symmetric 
2-D Gaussian convolution kernel with standard deviation a. 

^ which indeed make a significant difference in FBP performance. 
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Fig. 14: Reconstruction in a ROI. Upper left: true image, 
upper right: AFBP on truncated projections (SNR = 22.05 
dB). Lower left: FBP on truncated projections (SNR = 15.73 
dB). Lower right: FBP on full data (SNR = 22.28 dB). 

Given a training set Ttr = {/i, /2, /at}, we define 
1 ^ 

Cx = ar^rnin = ^ ^ ||XR/, - * fjh (16) 

In words, the blur measure Cx is the standard deviation of 
a Gaussian kernel which action on images of the training set 
resembles the most (in MSE sense) the action of XR. Notice 
that the value of Cx is, up to a constant, the Full- Width- 
Half-Maximum (FWHM) of the effective Gaussian blurring 
kernel, since the FWHM of a Gaussian is 2.35 multiplied by 
its standard deviation. 

Notice that for the AFBP transform T trained via equation 
([T5]) with corresponding standard deviation ctt, we have Ct = 
cfT' For any other transform X the blur measure can be easily 
approximated by evaluating the optimizer of ([T6l) via a grid 
search on values of a. 

The blur measure will allow us to parameterize the se- 
quences {T^},{Si} of reconstruction algorithms in order to 
correctly compare their performance: for a given level of 
noise present in the sinogram, the output quality is a function 
of the blur measure applied in the reconstruction process. 
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Fig. 15: Sinogram completion method used for FBP. The hori- 
zontal band in the middle of the image represent the available 
projections, truncated to the central disc in image domain. 
Other rows are replica of the upper and lower margins of the 
band. This simple step helps improving the FBP performance 
substantially. 



More importantly, this notion allows us to unambiguously 
characterize the sequence of reconstructors, employed in the 
first stage of the SPADES algorithm, instead of specifying the 
parameters of apodization window for FBP or the values of 
(Ji for AFBP 

In the following numerical experiment, a sequence {T^j^ 
was trained on 15 geometric phantoms by optimizing the 
equatior[T5l in the setup described in the beginning of the Sec- 
tion |^Dl Using these phantoms, the values Ctj are computed 
empirically, as well as the values of Cs^ for FBP transforms 
with corresponding cut-off frequency qj (see equation ([H])). 
We used different values of the parameter p, to achieve best 
available quality for the FBP. All algorithms are applied to 
noisy truncated projections, mended by the aforementioned 
completion technique (see Figure [15]). In Figure [161 we present 
graphs of SNR values of phantoms reconstructed from noisy 
and truncated projections, with the same setup as in the 
previous experiment. Specifically, we plot average SNR values 
over a set of phantoms, as function of the blur measure. The 
comparison is carried out between the AFBP and the apodized 
FBP with different values of p. It can be seen that with 
p = 0.5, FBP achieves maximal quality at C ^ 2.1, which is 
substantially lower than SNR values achieved by the trained 
kernel of AFBP at C 1.2. 




0.5 1 1.5 2 2.5 3 3.5 

Blur measure ^ 



Fig. 16: Graphs of reconstruction quality versus the blur 
measure Ct, computed for the sequence of AFBP transforms 
and also for three sequences of apodized FBP transforms with 
different degrees p of Butterworth window. 



VI. SPADES IN THE ROI 

A. The Algorithm 

We return to the SPADES algorithm (see Section HVl) with 
AFBP transforms. Now it can be extended to the setup of ROI 
reconstruction from truncated projections. To that end, AFBP 
transforms of two types, corresponding to objective functions 
([T4|),([T5]) are used: 

(1) One transform T^, trained with ([T4l) . will produce an image 
of best quality (in SNR terms) attainable with AFBP at the 
given noise level, for images similar to those of the training set. 

(2) A sequence of transforms {T^K^^ trained via the equation 
([T5l) , will produce a set of image versions with a growing blur, 
corresponding to standard deviation values {(7^}^^^}. we use 
/ = 10 such transforms. 
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The training of AFBP versions is performed on a given 
set of images. Afterwards, the same training set is used to 
compute the weights of the neural network, while the inputs 
are obtained from the linear AFBP image estimates. When 
all components of the SPADES reconstruction are ready, the 
processing of the sinogram data Qn for an unknown image is 
carried out in three steps (see a diagram in the Figure [5]): 

1) Compute preliminary images fi = Tignji = 1,2,...,/, 
and fo = Togn- 

2) For each pixel p in the ROI, build a vector Xp of features 
from the corresponding values of fi , fo (see Section 
HV-Al for details). 

3) Use the neural network to compute the value of the 
output image at p from the vector of features Xp. 

B. Computational Complexity of SPADES 

The computational complexity is calculated for images of 
size nxn and for / + 1 preliminary reconstruction operators. 

(1) The standard FBP algorithm, which consists of a 1- 
D filtering of the sinogram columns and a Back-Projection 
transform, requires O(n^) computations for each of these 
steps. The AFBP differs in the 2 — D sinogram filter, which 
uses about 3 — 5 times more operations than sl 1 — D filter, 
and in a 2-D post-processing filter which consumes O(n^) 
computations: we used a 2 — I) kernel of dimensions ^/nx^/n, 
therefore its action requires n operations per pixel. Therefore 
AFBP has the overall complexity of O(n^) computations. 

(2) SPADES involves / + 1 versions of the AFBP transform. 
If the computing machinery allows to execute the / + 1 AFBP 
algorithms in parallel, time consumption of the preliminary 
stage equals to a single FBP reconstruction. Otherwise, the 
preliminary step requires 0{I ■ n^) operations. Typical value 
of /, used in our experiments, is / = 10. 

(3) The neural network involves {I -\-ll)N weights, which 
participate in the computation of every pixel in the image. 
Here TV is a number of neurons and (again, in our experieece) 
it should be roughly equal to the size of the vector of features, 
/+10. Thus, computation of each intensity value for the output 
image requires ~ (/ + 10)^ ^ 400 operations. 

(4) If the number of pixels in ROI is M, the overall 
complexity is O((/+l)n^ + (/+10)2M) - O(10n^ + 400M) 
operations. 

C Numerical Experiments - SPADES on the ROI data 

In the first numerical experiment, we used a training set 
of 15 geometric phantoms, described earlier. The chosen 
blur measure values Ct^ are equally spaced in the interval 
[0, 3.5]. Linear reconstruction of the ROI in a test image with 
AFBP algorithms was conducted from noisy and truncated 
projections (see experiment setup in Section IV-DI The neural 
network was applied on the resulting images to produce the 
final output. 

Outcomes of three different algorithms are visually com- 
pared in Figure [T71 a linear reconstruction scheme AFBP, a 
locally-adaptive direct algorithm SPADES and a statistically 
based iterative algorithm [6] (see the Appendix section for de- 
tails of our implementation). It can be observed that SPADES 



has succeeded in reducing the background noise, stemming 
from low X-ray intensity, restoring the piece-wise constant 
texture of the phantom. The SNR value of SPADES output 
is higher by 1.15 dB than the SNR of best linear algorithm 
output we could produce. The statistical reconstruction result 
has sharper edges and more homogeneous ellipses, at the 
price of higher computational time. Notice that because of 
artifacts its SNR value is still lower than that of SPADES 
output. The reconstruction was further carried out on 23 other 
test images. The average SNR values of the ROI quality are: 
FBP (on truncated projections) - 15.80 dB, FBP (on full 
data) - 18.89 dB, AFBP (on truncated projections) - 19.31 
dB, SPADES (on truncated projections) - 20.30 dB. Such 
quality improvement was possible due to the fact FBP displays 
discretization artifacts, which restrict the reconstruction quality 
even in the absence of noise and with no truncation of the 
projections. The learned kernels of AFBP succeed to partially 
compensate this drawback. 




Fig. 17: Upper left: true image. Upper right: best AFBP 
version (18.44 dB). Lower left: SPADES reconstruction (19.81 
dB). Lower left: Statistical reconstruction (19 dB). 

The experiment was repeated for a ROI reconstruction in 
the Forbild head phantom, designed in the framework of the 
Forbild project 0. A set of ten similar phantoms, obtained from 
the head phantom by random perturbations of the radii and 
locations of the composing ellipses, served as a training set 
for SPADES. The behavior of the involved algorithms (AFBP, 
SPADES and the statistical reconstruction) was similar to one 
observed in the Figure [T71 Numerically, the quality of ROI 
recovered by AFBP in the Forbild phantom was 17.93 dB, 
statistical reconstruction has achieved the same value, and 
SPADES has demonstrated 20.1 dB. Due to the lack of space, 
we omit the accompanying graphics. 

Another experiment was carried out with the clinical CT 
images, described earlier in the Section IV-DI The AFBP 
output, generated in the earlier experiment on CT data (Section 
IV-DI Figure [141), is used, along with a sequence of blurred 

^http://wwwimp.uni-erlangen.de/forbild/english/results/head/head.html. 
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versions of the test image, as the input data for the neural 
network. The parameters of AFBP , as well as the neural 
network, were trained using a set of 10 clinical images (some 
are presented in Figure [12]). Output of AFBP, SPADES and 
the statistical algorithm are displayed in Figure [TSl Their 
behavior is similar to the one observed on the phantoms. 
The linear AFBP produced a goo(£| preliminary image, which 
nevertheless contains noise propagated from the projections. 
SPADES succeeds to remove the noise substantially, resulting 
in an image close to the ground truth. The statistical algorithm 
also produces a clean image, but the artifacts introduced by 
the Huber prior deteriorate its quality. SNR values reflect on 
these observations. 

In the reconstruction of 16 additional test images, average 
SNR gain of the SPADES algorithm over the AFBP was 
1.93 dB. This shows consistency in the behavior of the two 
algorithms, visually displayed in the Figure [TSl 




Fig. 18: Upper left: true image. Upper right: best linear AFBP 
(22.04 dB). Lower left: SPADES reconstruction (24.58 dB). 
Lower left: Statistical reconstruction (19.24 dB) 



VIL Discussion and Conclusions 

We have proposed and demonstrated a CT reconstruction 
algorithm based on a local fusion of linear estimates of the 
sought image. The linear reconstructors were designed to 
automatically adapt to the specific family of images and the 
noise level, as well as the missing data setup. Resulting direct 
non-linear reconstruction algorithm produces images inher- 
ently different from those obtained with EBP or other linear 
methods: the noise level and artifacts are substantially reduced, 
while the edges present in the image are well preserved. We 
notice that the proposed method can in principle improve any 
given reconstruction method, by fusing its output image with 
a number of other image estimates by an adequately trained 
neural network. 

The proposed scheme involves very few parameters, mostly 
the design details - number of linear reconstructors and their 

^see section IV-D] for comparison with FBP 



blur measure values, structure of the vector of features and 
number of neurons. All these are easily tuned in a practical 
setup of specific implementation of the algorithm. Admittedly, 
the neural network requires a delicate treatment, but in our 
experience (and especially for specialists in Machine Learning) 
it is fairly manageable. In fact, we used the simplest one-level 
feed-forward model of the neural network, and it is plausible 
that using more complex network (with larger training set) will 
improve the algorithm behavior. 

We notice that the proposed scheme performs an isotropic 
processing in the image domain: the inputs of the neural net- 
work are preliminary images fi = TR/ which approximate 
the true image f{x) convolved with a radially symmetric 
kernel. The FBP algorithm also has a radially symmetric 
PSE. A promising future direction of work is an extension of 
SPADES to an anisotropic algorithm, where elliptical kernels 
with varying directions and dimensions are used. The neural 
network should automatically choose the relevant direction 
(or combination of some), to use the appropriate preliminary 
image version at each point 

All the numerical results presented in this paper were 
obtained from experiments executed entirely in the Matlab 
environment. The "ground truth" images are projected using 
one specific implementation of the discrete Radon transform, 
and also used as the ideal reference for the training pur- 
poses. Admittedly, the results are affected by this setup, 
hence the experiments should be rendered as synthetic. The 
proposed methods can be extended to a practical setup by 
using mechanical phantoms with geometric details and their 
discrete computerized models. Despite the absence of realistic 
clinical data experiments and industrial implementations of 
the FBP in the experiments, the differences observed between 
the existing and the proposed techniques are consistent and 
explicit. Therefore, we hope that presented techniques will 
also exhibit similar behavior when applied in the industrial 
CT setup. 

Notice that incorporating the proposed algorithm into an 
existing CT machinery will only require minor software modi- 
fications - replacement of FBP filters and addition of a learning 
machine. Since SPADES is not slower than the standard FBP 
algorithm (when parallel computation of FBP instances is 
available), it represents an appealing alternative for full- scan 
reconstruction and especially for ROI recovery from partial 
data. 

In a broader view of numerical algorithms for inverse prob- 
lems, the proposed technique has quite a general underlying 
principle: a fusion of a sequence of solutions with a certain 
structure, using a discriminatively trained learning machine. 
One possible application of this principle can be considered 
for the very common method of solving an inverse problem 
by optimization of the objective function of form 

X* = arg mm{D{x, y) + \P{x)}. (17) 

X 

Here is a set of measurements, x is the sought signal, D{x^y) 
is the Data component which encourages the consistency of the 
solution with the observations and P{x) is a Prior component, 
expressing the beliefs on the nature of the signal x. The 
optimizer is usually obtained by an iterative update of the 
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initial solution xq. The parameters of the algorithm, such as 
the value of A and the number of iterations, are crucial for 
successful reconstruction but are difficult to estimate and often 
are data-dependent. 

An appealing solution to the problem of chosing the parame- 
ters would be to compute a set of optimizers x* , corresponding 
to different values of A and/or the number of iterations. These 
preliminary optimizers would then be fused into the final out- 
put by means of a neural network (or other learning machine), 
trained on samples from a training set. Such approach (on 
expense of higher computational complexity) is expected to 
produce a solution x of quality higher than any individual 
X*, as was the case in numerical algorithms presented in this 
paper. 

Appendix 

A. Generating the Poisson Noise 

The Poisson noise is introduced into photon counts yi as 
a consequence of limited values of yi themselves; the lower 
the total count of the photons, received by a detector per unit 
of time, the smaller the accuracy of the measurement for the 
corresponding line integral. The mathematical model of such 
noise is the Poisson distribution of photon counts: yi is an 
instance of the random variable F^, 



and Rs{f) is an edge-preserving penalty which encourages 
image smoothness: 



Yi ~ Poisson{Xi)^ A^ = Iqc' 



(18) 



Here /o is the scanned image (a 2-D slice of an object), Iq is 
the photon count in absence of obstacles (proportional to the 
initial X-ray intensity), and (R/)i is the line integral of the 
image along the line corresponding to detector i. 

In our synthetic experiments, we first compute R/ and 
then generate the set of photon counts. This is done as 
follows: (1) Denote g = R/ and compute y^ = /oe~^^^^ for 
every bin i in the sinogram. Naturally, the maximal photon 
count is max^{?/^} = Iq since min{g) = 0. The minimal 
count ymin = ^^^iiVi} depends on the scaling of g. (2) 
We chose this scaling such that ymin = 60, in order to 
avoid the problematic numerical behavior related to a low 
photon count, contaminated with the Poisson noise. Notice 
that the ratio lo/ymin should be constant when different 
values of Iq are used. (3) Draw instances of random vari- 
ables yf ~ Poisson{y^), index i runs over all the bins in 
the sinogram. (4) convert back to the sinogram domain, by 
computing gn = -log{y'^ / Iq). 

B. Statistically -based Iterative Reconstruction 

We implement a Penalized-Likelihood reconstruction algo- 
rithm, defined in [6] for mono-energetic X-ray scan. Given the 
noisy measurements y of photon counts, the output image / 
is computed by minimizing the Penalized Likelihood equation 

^{f) = L{f\y) + m6{f) (19) 

where L{f\y) is the negative log-likelihood of the Poisson 
distribution. 



N 



L{f\y) = ^{^i- ydog{\i) + log{yi\)} , A^ 



(20) 



(21) 



with D/ = [Va^/, Vyf] consisting of the directional derivative 
maps of the image /. We compute the derivatives by the 
central difference approximation. It is defined for x-axis by 
fx{px,py) = f{px + l,py) - f{px - l,py) and similarly for 
the ?/-axis. The Huber penalty ^1)5 is defined by 



x^ 

T 

8\x\ 



2 



if \x\ < 8, 
if \x\ > 6. 



The minimizer of equation ([T9b is computed using a Gra- 
dient Descent with a simple line- search function (we did 
not pursue a computationally effective implementation, since 
only the quality of the algorithm output is considered). The 
parameters /3, S and the number of iterations were tuned on 
the training images for each image to attain the best possible 
SNR value. 
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